Continuous Marker Prognostic Power

Description

This program uses simulation to approximate the power of a test of a marker association with survival data assuming an exponential regression model. The calculations include a proportional hazards test using the model where the log-hazard ratio is linear in a continuous marker distribution. The choice of distributions include uniform, normal, and lognormal. All covariates are mean centered for the analysis. The magnitude of the effect size can be specified based on a unit of the prognostic factor or the interquartile range of the marker.

It is useful to note that if the underlying distribution is assumed to be uniform, and the proportional hazards model regression log-hazard ratio is linear in the marker, then effect size associated with the interquartile range will correspond to approximately same effect as a marker dichotomized at the median (or other quantile).

In addition to power calculations, the approximate confidence intervals of parameter effects and expected number of events are also presented.

The impact of the continuous prognostic model assumptions in terms of the expected events can be helpful in gaging assumptions on expected events relative to the overall trial design.

Input Items

The user is prompted to enter values for the following items. For items that have initial default values set, the default values are given in parentheses.

• Number of Simulations (1000): The number of simulations to run. For better approximations, increase the number of simulations.
• Total Sample Size
• Alpha Level (0.05): \(\alpha\), the 1 or 2-sided significance level
• Accrual Period: Length of the accrual period
• Follow-up Period: Length of the follow-up period, i.e. the time from end of accrual to analysis
• Hazard Rate (marker value = 0): The hazard rate at marker value equal to 0
• Effect Size (Hazard Ratio per Unit): The hazard ratio associated with a one unit change in the marker
• Unit (Interquartile Range): Unit of marker effect is either Interquartile Range or Unstandardized
• Sided: One-sided or two-sided test
• Marker distribution (Uniform): Either standard uniform, or normal (standardized to a mean of zero), or lognormal
• Log-Normal Distribution Parameter (1): The spread/shape parameter (sigma) for each of the distributions. It is fixed to 1 for uniform and normal and defaulted to 1 for log-normal. The log-normal density is specified as \[ f(x) = \frac{e^{-(log(x) - \mu)^2/(2\sigma^2)}}{\sqrt{2\pi}{\sigma}x} \] where μ and σ are the mean and standard deviation of the logarithm of the normal random variable.

Output Items

• Estimated power
• Estimated effect and empirical histogram simulated expected confidence interval
• Expected number of events in the entire study

Note

This is a simulation-based calculation and therefore may take several seconds to calculate depending both on the number of simulations and sample size.

Statistical Code

The program is written in R.

View Code

function(accrual, follow, n, alpha, sided, lambda1, unit, HR1, dist, dist_param, nsim) {
    # Power for continous marker #
    # Note that covariates are mean centered in calculations #

    library(survival)

    beta1 = log(HR1)

    # This is where the biomarker distibution is chosen
    if (dist == "unif") { dist.generate = runif }
    if (dist == "norm") { dist.generate = rnorm }
    if (dist == "lnorm") { dist.generate = rlnorm }

    testscorea = coefa = nevent = rep(NA, nsim)

    # Loop over simulations 
    for (j in 1:nsim) {

        # Generate Data #
        #---------------#

        biomarker = dist.generate(n, 0, dist_param)

        biomarker0 = biomarker - mean(biomarker)
        if (unit == "IQR") {
            biomarker0 = biomarker0 / IQR(biomarker0)
        }

        tt = (-log(runif(n))) / exp(log(lambda1) - beta1 * biomarker0) # generate true survival by taking exponential(1) / exp(alpha + beta * biomarker)
        cc = follow + runif(n) * accrual                               # generates true censoring times
        status = 1 * (tt < cc)                                         # generates 0,1 censoring indicator
        times = tt * status + cc * (1 - status)                        # observed failure times (times under observation)

        # Do Analyses #
        #-------------#

        rega = coxph(Surv(times, status) ~ biomarker0)                                 # do cox regression continuous
        coefa[j] = rega$coef                                                           # cox reg coefficient (for continuous variable)
        testscorea[j] = rega$score                                                     # score test
        nevent[j] = sum(status)
    }

    powercont = mean(testscorea > qchisq(1 - alpha * sided, df=1))
    coefcont = quantile(coefa, c(alpha / sided, .5, 1 - alpha / sided))
    events = mean(nevent)

    if (dist == "lnorm") {
        plot_name = "log_normal_distribution.png"
        png(filename = plot_name)
            x1 = qlnorm((1:50) / 50, meanlog = 0, sdlog = dist_param)
            y1 = dlnorm(x1, meanlog = 0, sdlog = dist_param, log = FALSE)
            plot(x1, y1, xlab = "x", ylab = "log normal density", type = "n")
            lines(x1, y1, lwd = 2)
            title("Log Normal Distribution") 
        dev.off()
    } else {
        plot_name = ""
    }

    result =
        list(power_continuous = round(powercont, digits = 3),
             expected_events = (floor(events * 10) / 10),
             effect = round(exp(-coefcont[2]), digits = 3),
             lower_ci = round(exp(-coefcont[3]), digits = 3),
             upper_ci = round(exp(-coefcont[1]), digits = 3),
             session_key = strsplit(getwd(), "/")[[1]][4],
             log_normal_plot_name = plot_name)
    return(jsonlite::toJSON(result, pretty = TRUE))
}